Methods

library(tidyverse)
library(ggdist)
library(ggside)
library(easystats)
library(patchwork)
library(brms)
source("functions.R")

df <- read.csv("../data/study2.csv") |>
  mutate(
    Date = lubridate::dmy(Date),
    Participant = fct_reorder(Participant, Date),
    Screen_Refresh = as.character(Screen_Refresh),
    Illusion_Side = as.factor(Illusion_Side),
    # Illusion_Effect = fct_relevel(as.factor(Illusion_Effect), "Incongruent", "Null", "Congruent"),
    Illusion_Effect = fct_relevel(as.factor(Illusion_Effect), "Incongruent", "Congruent"),
    Block = as.factor(Block),
    Education = fct_relevel(Education, "Doctorate", "Master", "Bachelor", "High School", "Other")
  )

# Fix precision
for(ill in unique(df$Illusion_Type)) {
  data <- df[df$Illusion_Type == ill, ] 
  i <- 10
  while (length(sort(unique(signif(data$Illusion_Difference, i)))) != 8) {
    i <- i - 1
  }
  df[df$Illusion_Type == ill, "Illusion_Difference"] <- signif(df[df$Illusion_Type == ill, "Illusion_Difference"], i)
  if (i != 10) {
    message(ill, ": Illusion_Difference values != 8. Rounded to ", i, ".")
  }
}

# Transformation
df <- df |> 
  mutate(    
    Illusion_Difference_log = log(1 + Illusion_Difference),
    Illusion_Difference_sqrt = sqrt(Illusion_Difference),
    Illusion_Difference_cbrt = Illusion_Difference**(1 / 3),
    Illusion_Strength_log = sign(Illusion_Strength) * log(1 + abs(Illusion_Strength)),
    Illusion_Strength_sqrt = sign(Illusion_Strength) * sqrt(abs(Illusion_Strength)),
    Illusion_Strength_cbrt = sign(Illusion_Strength) * (abs(Illusion_Strength)**(1 / 3))
    )

Exclusions

# plot(estimate_density(dfraw[dfraw$Participant == "60684f29dbfe1bb2059e5e27_rkqoy", "RT"]))

# Dear participant, thank you for participating in our studfy. Unfortunately, we didn't receive your data :( did something happen? Some technical issue? We would like to kindly ask you to return your participation so that we can open up more slots. Thank you in advance, and apologies for the inconvenience! 

# Dear participant, thank you for participating in our study. Unfortunately, our system detected multiple issues in your data (such as implausibly short responses, and random-like pattern of answers), which makes it unusable. We understand that you might have been in a hurry or had some other issues; we hope to open-up more slots in the future would you be interested to participate again.

# Dear participant, thank you for participating in our study. Unfortunately, our system detected multiple issues in your data (such as implausibly short responses, and random-like pattern of answers - in particular in the 2nd part of the study), which makes it unusable. We understand that you might have been in a hurry or had some other issues, and so we kindly ask you to return your participation; we hope to open-up more slots in the future would you be interested to participate again.

# Just received the results: in your case, the three most prominent issues were your response pattern that was equivalent to random (in the sequence of keystrokes) that led to a chance-level performance (that was also significantly different from the rest of the population). Moreover, your reaction time distribution was also very different from the norm, with a vast majority of implausibly short responses (i.e., that are faster than the time it takes the brain to process a visual input). This issue was even more prominent in the second block (after the break), which typically happens when participants are in a rush to finish. Finally, your overall completion time was also significantly below the average. Again, we apologize, we understand that your time is valuable, but unfortunately we run too on limited funds and can hardly spare more on unusable data. Sorry for the inconvenience, we will try to open-up more slots soon.

outliers <- c(
  # Error rate of 48.8% Very short RT
  # Prolific Status: REJECTED (06/08)
  "5e66ed53de7896486f9a71f8_p8ckk",
  # 2nd block of responses very fast
  # Prolific Status: REJECTED (15/08)
  "61443e47c248a375c899a0cf_9qz3u"
)

partial_outliers <- c(
  # 2nd block a bit bad
  "5c43cd414fe4f800016e4983_0zp37", 
  # Entire 2nd block bad
  "61564e378e974bdbe42763a2_hhufm")

We removed 2 participants upon inspection of the average error rage (when close to 50%, suggesting random answers) and/or when the reaction time distribution was implausibly fast.

For each block, we computed the error rate and, if more than 50%, we discarded the whole block (as it likely indicates that instructions got mixed up, for instance participants were selecting the smaller instead of the bigger circle).

Error Rate

dfsub <- df |>
  group_by(Participant) |>
  summarize(
    # n = n(),
    Error = sum(Error) / n(),
    RT_Mean = mean(RT),
    RT_SD = sd(RT),
    IPIP6_SD = mean(IPIP6_SD),
    PID5_SD = mean(PID5_SD),
  ) |>
  ungroup() |>
  arrange(desc(Error))

data.frame(Participant = c("Total"), t(sapply(dfsub[2:ncol(dfsub)], mean, na.rm=TRUE))) |> 
  rbind(dfsub) |> 
  knitr::kable() |>
  kableExtra::row_spec(1, italic = TRUE) |> 
  kableExtra::row_spec(which(dfsub$Participant %in% outliers) + 1, background = "#EF9A9A") |> 
  kableExtra::row_spec(which(dfsub$Participant %in% partial_outliers) + 1, background = "#FFECB3")

Participant	Error	RT_Mean	RT_SD	IPIP6_SD	PID5_SD
Total	0.220	793	601	20.34	20.48
5e66ed53de7896486f9a71f8_p8ckk	0.488	311	337	NA	NA
61443e47c248a375c899a0cf_9qz3u	0.359	417	434	19.93	21.55
60d60bb9ced0566454d42750_020cx	0.349	618	227	22.37	25.55
5c43cd414fe4f800016e4983_0zp37	0.339	796	657	23.26	12.29
60d84da28f38998c15e2a4fa_xw2cl	0.333	702	376	25.12	31.38
5cbb77bf418c540001641e67_sp53j	0.324	870	547	23.39	15.09
60c5ae1d153847577f543609_u9w9l	0.302	878	950	32.58	46.24
611d69ffd22b95ca500604b2_hs0kt	0.300	855	528	24.03	15.30
6140b5b201e55201433c2383_60sey	0.299	674	303	9.79	16.05
5f7cd5bc00f2a016f22e74d2_4qe89	0.296	521	269	25.43	36.14
61564e378e974bdbe42763a2_hhufm	0.289	663	645	25.81	25.94
5c067e6329abc0000199ae45_lbvy6	0.283	579	178	14.20	17.66
60f1ce9c530c91c0a6e9dab0_1rwja	0.283	766	399	22.55	23.73
5efbc3a84e1b6416ad27a3ec_qdax2	0.281	837	2198	14.87	11.78
60e72879166656b7391b2eb2_4t9nj	0.280	897	512	19.71	25.95
6059025fd396ddcd93fe50aa_pytn9	0.279	1112	1363	20.25	23.23
613bb8d7741ae6ca04e0431b_gk94a	0.277	921	565	NA	NA
5cb632bab9c4860001700fe6_dctvu	0.277	719	570	16.66	14.08
613d8489ab93bd9635c12051_zwbj3	0.272	730	269	31.39	23.54
60560efe16c645d2a9bd0daa_llvby	0.270	744	301	22.17	18.34
613c1ac80e6fc9ac43eac753_dbdal	0.266	743	396	20.36	18.70
6170adf5330be4e266555619_68ago	0.262	824	542	18.08	13.53
614b55e22ff3944a165736bb_v16g9	0.261	518	170	11.49	14.03
5eaef8702b68455d6e130595_wh7q8	0.259	707	338	14.92	12.77
606781fa7584b511b845d52c_cu5nv	0.256	550	152	NA	NA
5f9aba6600cdf11f1c9b915c_at42y	0.254	788	593	22.87	20.80
5f2e1510b83f09172c62ba16_22y38	0.251	1023	1708	8.40	16.45
60feee53991033aa986fa755_lp5dv	0.251	658	330	13.83	17.35
5bdbfff9ee652a0001efca64_vvdyu	0.248	627	302	23.10	19.36
61438c83c0d9428ad7b4e4b8_t0wjb	0.248	731	336	22.27	34.83
612383021976c4e2cac79581_0fzdd	0.247	738	235	30.22	23.96
60087e9189ddb34b79b08be5_pd14q	0.245	613	202	26.64	18.85
589630d22a697d0001cfca12_buz00	0.243	850	253	18.60	30.58
5f46a8269f9d230534df642a_gql1v	0.242	1354	1254	12.23	10.61
5c437f6a4fe4f800016e3d52_a73wb	0.241	554	177	30.65	25.55
611fb1635c5f446a6bb2848d_6fcvf	0.241	990	513	37.35	28.90
5f494676e0837811e94a4007_qfzu1	0.240	940	670	16.18	21.86
5f983f3b00bbab0c4b28149a_l7hnq	0.240	646	225	11.55	21.11
5d3c6e745602310001bca8aa_6bdtb	0.238	662	495	NA	NA
61645edbf8a9840feeb735b6_qot1q	0.238	1410	1660	17.48	16.72
614060a52d7c64c27ef9887c_43ppx	0.238	794	628	18.33	13.04
611e6b43e936390a65656ecd_rr1n5	0.238	578	166	27.05	19.75
5b8c8e7d0c740e00019d55c3_ejsw9	0.236	674	313	15.44	17.20
5e98b398e4a84f0eb3755e51_yy8z9	0.235	574	241	13.95	23.42
5c8fcb13dd3b360015b4dff6_dvbq9	0.234	643	327	NA	NA
5f6b027e6eae971e2fa13594_5qh7m	0.233	555	180	16.80	13.21
5f304d07956a3f48639e3aa6_03qfs	0.231	698	286	26.22	27.93
615d5311a9fbf1f98aa62a5a_n2x53	0.229	495	182	18.67	20.43
610cda1332fb63830158c55d_95zqk	0.228	691	458	24.55	24.63
60ec51c51a3158a50ded8a3e_jj888	0.227	676	252	14.28	10.15
613bab4160cd2fe41f9b1a27_kepjv	0.226	587	292	20.87	22.09
60efc61b7b9f908aea94ed57_z9d93	0.224	1506	12083	26.40	17.62
612cb8236aa1cef2599e6f59_6jam2	0.221	622	751	18.93	16.50
60d1224bb86fedf1819ef67b_tvh43	0.221	790	579	17.25	26.71
6166c19fb69f46ba9749489b_6ym17	0.220	735	313	22.95	25.24
60142d79e4c5e23b4730a69e_98qsp	0.220	618	189	17.14	19.61
606f598baf5e0ffda45e9e47_qayzo	0.220	1327	1557	9.45	14.82
611ee495b12c10171def6816_9oz56	0.219	1226	928	26.97	25.54
612819d5c4f57a5fb23d03b3_sc9zu	0.218	630	467	14.06	16.49
611183c67911657721712872_sdh4z	0.218	846	529	14.36	12.28
60fa8b48cbeec7cf62317c21_glept	0.216	718	346	27.40	36.23
5c1bb460a05a64000125c522_ohxap	0.216	564	169	30.73	20.37
60feabd18109e08e594540f8_r9m9e	0.216	622	191	23.11	29.01
5e8e06cbdf10df0a811b6f29_g2tm6	0.215	610	293	NA	NA
6112810f76d46087c3e2d208_51bwf	0.215	770	663	22.42	20.99
5e283bdff5bb81010fcc6cde_9zgbc	0.214	1446	1197	13.93	22.34
5d2508273854c300016c04aa_1qm3y	0.214	757	353	11.82	13.46
5b232f6838fc0c000131438c_g758c	0.212	1024	766	28.16	21.16
5e938e072a4ce9077f8157b9_egf3b	0.212	1012	781	16.75	27.39
5c42329384ef0f00011882ca_tae13	0.209	633	268	27.45	27.36
60c0fe89173c2b1f4fe2c578_a23gw	0.209	783	341	31.55	30.88
60f5fa488f7b1381d175ecb5_qwvfs	0.208	670	341	20.86	14.25
6010b797bbe6be440425665c_8qarc	0.207	816	521	11.95	15.38
605d5ffca1324379006a3599_703q5	0.207	512	139	11.48	11.99
5e6e14f53c76d23b934a67f3_ghxwv	0.206	911	2764	33.56	27.03
6139133041546096188da709_kr98d	0.206	799	465	21.71	19.13
5fcbd87e70a25c99595f012d_dkj0n	0.205	564	155	17.76	20.20
61602320c1aec7e6f8b5e3f4_adnrd	0.205	558	177	22.12	8.45
60c21f5f1389f65d2d88a3bf_lxzqo	0.204	1030	466	35.66	20.25
5f264332488e13372281c1c8_b082g	0.204	745	359	19.00	11.51
60a5404314ae1a9d5b88565d_fd4x8	0.204	763	274	15.46	15.63
61264d8a5de54eb1e42ed381_6x9p2	0.203	735	242	25.56	6.79
60009275300e1c14bf67ce83_7n71t	0.203	602	202	22.41	31.03
5fa18c4a48882016b4b143f7_kn4c0	0.202	965	984	20.97	32.36
62d6d791bb6448ae52929ebd_ncqkx	0.202	776	467	30.03	27.75
61395b705a38ede105ed7982_jkfrp	0.202	639	298	20.92	25.94
5f439ac86b921f5fac4f55ad_afpfo	0.196	884	656	22.10	16.21
60127f78468af02bcca0b7a0_lad8d	0.196	576	250	NA	NA
5f97e6601f6d0e016087fc91_bac1g	0.195	649	223	20.92	19.93
614f8321333460db43b7f229_mpxhl	0.194	1645	983	NA	NA
5f733312e59fd2000a74b8cf_lbgw8	0.194	972	612	12.46	16.09
5fc78656c7368206ab1c5174_4t6uk	0.192	640	306	13.37	22.43
614b23f447bd48a3801fbf95_dlqf9	0.190	739	437	16.14	18.70
6106ac34408681f3b0d07396_y2n16	0.189	613	221	27.84	24.74
60e866921fd3ba5f3007c2aa_7toga	0.189	987	1084	16.00	23.38
5ed7a7a467a98224295459ff_0l6hx	0.187	531	107	13.22	10.44
5ceef0cd45c1ec001920d548_kurpx	0.185	1062	758	22.10	22.66
610af9910c4d079d84e855f4_z8j00	0.185	764	360	25.17	16.11
5e6d381c029fa42f9f48c105_kcgfo	0.185	869	562	24.09	20.61
5f782504951dda38db3ef0b6_qsfgg	0.184	1072	798	14.89	15.58
613646a48d30fb040cf5ffed_g5q2e	0.184	702	184	14.16	15.04
614a4c9d8172751921452626_mc4l4	0.184	772	305	31.56	28.66
5fdfd04b9bf07d83b2e5f780_qcv0k	0.183	886	537	23.91	23.77
613231ee733769a9aff9c6c5_ps1tt	0.180	896	619	11.43	25.06
60d26bfdacc7efc851571e74_l8k5r	0.180	718	372	18.15	27.52
5eb3a734d249ac18a413063a_ec0gy	0.179	921	527	20.16	14.07
5f61341fe5033f0b6327fb7f_lg9vw	0.179	775	384	13.89	20.64
5ede636d21cf560bf7aeb29f_03zzd	0.178	741	532	17.33	21.62
6026f1c19366fc4778bfd055_nuccw	0.173	873	568	15.28	19.65
60ddfb3db6a71ad9ba75e387_9elmo	0.173	692	287	16.89	12.46
60d340fe43207d6a5bc89ca1_m32vm	0.173	1090	1057	14.99	18.71
5f31de301a469e11624f7948_3n2lv	0.164	783	311	19.53	26.15
611bb2a6d9678117db737236_pf2qy	0.163	904	549	12.64	12.23
6147c21c933739516362d0ee_wcaat	0.163	900	460	24.99	24.26
5ecd58b69c382e0a177e0f61_p3yev	0.162	652	279	NA	NA
6067f370472020b9dc5cf71b_f2or0	0.162	861	337	16.49	17.12
61034f24da19cc56177b8b59_43rrv	0.161	805	418	28.19	13.51
614729093372b9e84bc48f8e_v4mcv	0.159	952	530	26.37	17.89
62da72f3fb2ff174670bc85f_vb5a9	0.157	786	365	18.68	13.44
5e8348450e128b0539d77ddc_z361c	0.154	820	523	9.39	24.93
5f8cbc5c355ea745e6cef2ca_5wo4f	0.152	765	406	10.78	23.51
610a9a0723dac03a63f15035_gclx4	0.151	795	312	28.52	35.59
612684c6d2cef3acadf5d77b_pv444	0.151	668	389	26.37	26.56
5d6e3252dfe6d10001dbe508_1aloo	0.150	816	379	25.32	19.75
60cefa69352cbf2549f2bf35_e13a8	0.150	688	221	24.42	25.38
5efd2964d36f63162f263795_xksu5	0.148	849	646	9.73	14.50
60ec6cb208c74c0d7e4794b9_hjay3	0.147	865	720	18.17	17.81
56cb8858edf8da000b6df354_m4821	0.146	579	211	12.57	10.43
6165e0b51a883e9db8cc7146_z0xax	0.145	816	647	23.25	14.69
5cdc88cc50f50f001a783112_tw4fr	0.118	1933	1570	23.15	13.20

Reaction Time Distribution

# RT distribution
p <- df |> 
  # filter(Date > "2022-08-12") |> 
  filter(RT < 10000) |> 
  estimate_density(select = "RT", at = c("Participant", "Block")) |>
  group_by(Participant) |>
  normalize(select = "y") |>
  ungroup() |>
  mutate(Participant = factor(Participant, levels=levels(df$Participant))) |> 
  mutate(color = ifelse(Participant %in% outliers, "red", ifelse(Participant %in% partial_outliers, "orange", "blue"))) |>
  ggplot(aes(x = x, y = y)) +
  geom_area(data = normalize(estimate_density(filter(df, RT < 10000), select = "RT"), select = "y"), alpha = 0.2) +
  geom_line(aes(color = color, group = interaction(Participant, Block), linetype = Block)) +
  geom_vline(xintercept = 125, linetype = "dashed", color = "red") +
  scale_color_manual(values = c("red" = "#F44336", "orange" = "#FF9800", "blue" = "blue"), guide = "none") +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  coord_cartesian(xlim = c(0, 3000)) +
  theme_modern() +
  theme(axis.text.y = element_blank()) +
  facet_wrap(~Participant) +
  labs(y = "", x = "Reaction Time (ms)")
# p
# ggsave("figures/outliers.png", p, width=20, height=15)

# Filter out
df <- filter(df, !Participant %in% outliers)

Error Rate per Illusion Block

temp <- df |>
  group_by(Participant, Illusion_Type, Block) |>
  summarize(ErrorRate_per_block = sum(Error) / n()) |>
  ungroup() |>
  arrange(desc(ErrorRate_per_block))

temp2 <- temp |>
  filter(ErrorRate_per_block >= 0.5) |>
  group_by(Illusion_Type, Block) |>
  summarize(n = n()) |>
  arrange(desc(n), Illusion_Type) |>
  ungroup() |>
  mutate(
    n_trials = cumsum(n * 56),
    p_trials = n_trials / nrow(df)
  )

# knitr::kable(temp2)

p1 <- temp |>
  estimate_density(at = c("Illusion_Type", "Block"), method="KernSmooth") |>
  ggplot(aes(x = x, y = y)) +
  geom_line(aes(color = Illusion_Type, linetype = Block)) +
  geom_vline(xintercept = 0.5, linetype = "dashed") +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  labs(y = "Distribution", x = "Error Rate") +
  theme_modern()

p2 <- temp2 |>
  mutate(Block = fct_rev(Block)) |>
  ggplot(aes(x = Illusion_Type, y = p_trials)) +
  geom_bar(stat = "identity", aes(fill = Block)) +
  scale_x_discrete(expand = c(0, 0)) +
  scale_y_continuous(labels = scales::percent, expand = c(0, 0)) +
  labs(y = "Percentage of Trials Removed", x = "Illusion Type") +
  theme_modern() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

p1 | p2



# Drop
df <- df |>
  group_by(Participant, Illusion_Type, Block) |>
  mutate(ErrorRate_per_block = sum(Error) / n()) |>
  ungroup() |>
  filter(ErrorRate_per_block < 0.5) |>
  select(-ErrorRate_per_block)

# Drop also participant with bad second block
df <- filter(
  df,
  !(Participant %in% partial_outliers & df$Block == 2))

rm(temp, temp2)

Reaction Time per Trial

df <- df |>
  group_by(Participant, Error) |>
  mutate(Outlier = ifelse(Error == 0 & (RT < 125 | standardize(RT) > 4), TRUE, FALSE)) |>
  ungroup()

p1 <- estimate_density(df, select = "RT", at = "Participant") |>
  group_by(Participant) |>
  normalize(select = "y") |>
  ungroup() |>
  merge(df |>
    filter(Error == 0) |>
    group_by(Participant) |>
    summarize(Outlier = mean(RT) + 4 * sd(RT))) |>
  mutate(Outlier = ifelse(x >= Outlier, TRUE, FALSE)) |>
  ggplot(aes(x = x, y = y)) +
  geom_area(data = normalize(estimate_density(filter(df, RT < 10000), select = "RT"), select = "y"), alpha = 0.2) +
  geom_line(aes(color = Participant, linetype = Outlier), alpha=0.2) +
  geom_vline(xintercept = c(125), linetype = "dashed", color = "red") +
  scale_color_material_d("rainbow", guide = "none") +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  guides(linetype = "none") +
  coord_cartesian(xlim = c(0, 4000)) +
  theme_modern() +
  theme(axis.text.y = element_blank()) +
  labs(y = "", x = "Reaction Time (ms)")


p2 <- df |>
  group_by(Participant) |>
  summarize(Outlier = sum(Outlier) / n()) |>
  mutate(Participant = fct_reorder(Participant, Outlier)) |>
  ggplot(aes(x = Participant, y = Outlier)) +
  geom_bar(stat = "identity", aes(fill = Participant)) +
  scale_fill_material_d("rainbow", guide = "none") +
  scale_x_discrete(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0), labels = scales::percent) +
  see::theme_modern() +
  theme(axis.text.x = element_blank()) +
  labs(y = "Percentage of outlier trials")

p1 / p2

We removed 1244 (0.78%) outlier trials (125 ms < RT < 4 SD above mean).

df <- filter(df, Outlier == FALSE)
df$RT <- df$RT / 1000  # Convert to second for better model convergence

Participants

dfsub <- df |>
  group_by(Participant) |>
  select(Participant, Age, Sex, Education, Nationality, Ethnicity, Duration, Break_Duration, Screen_Resolution, Screen_Refresh, Device_OS, starts_with("IPIP6"), starts_with("PID5")) |>
  slice(1) |>
  ungroup()

The final sample included 128 participants (Mean age = 26.1, SD = 7.1, range: [18, 69]; Sex: 46.9% females, 52.3% males, 0.8% other; Education: Doctorate, 0.78%; Master, 12.50%; Bachelor, 46.09%; High School, 34.38%; Other, 5.47%; Prefer not to Say, 0.78%).

plot_distribution <- function(dfsub, what = "Age", title = what, subtitle = "", fill = "orange") {
  dfsub |>
    ggplot(aes_string(x = what)) +
    geom_density(fill = fill) +
    geom_vline(xintercept = mean(dfsub[[what]]), color = "red", linetype = "dashed") +
    scale_x_continuous(expand = c(0, 0)) +
    scale_y_continuous(expand = c(0, 0)) +
    ggtitle(title, subtitle = subtitle) +
    theme_modern() +
    theme(
      plot.title = element_text(face = "bold", hjust = 0.5),
      plot.subtitle = element_text(face = "italic", hjust = 0.5),
      axis.title.x = element_blank(),
      axis.title.y = element_blank(),
      axis.text.y = element_blank()
    )
}

plot_waffle <- function(dfsub, what = "Nationality", title = what) {
  ggwaffle::waffle_iron(dfsub, what) |>
    # mutate(label = emojifont::fontawesome('fa-twitter')) |>
    ggplot(aes(x, y, fill = group)) +
    ggwaffle::geom_waffle() +
    # geom_point() +
    # geom_text(aes(label=label), family='fontawesome-webfont', size=4) +
    coord_equal() +
    ggtitle(title) +
    labs(fill = "") +
    theme_void() +
    theme(plot.title = element_text(face = "bold", hjust = 0.5))
}

p1 <- plot_distribution(dfsub, "Age", fill = "#FF9800")
p2 <- plot_distribution(dfsub, "Duration", title = "Total Duration", subtitle = "in minutes", fill = "#F44336")
p3 <- plot_distribution(dfsub, "Break_Duration", title = "Break Duration", subtitle = "in minutes", fill = "#3F51B5")

p4 <- plot_waffle(dfsub, "Sex") +
  scale_fill_manual(values = c("Male" = "#2196F3", "Female" = "#E91E63", "Other" = "#FF9800"))

p5 <- plot_waffle(dfsub, "Education") +
  scale_fill_viridis_d()

p6 <- plot_waffle(dfsub, "Nationality") +
  scale_fill_manual(values = c("South American" = "#FF5722", "Middle Eastern" = "#FF9800", "Western" = "#2196F3", "African" = "#4CAF50", "South African" = "#009688"))

p7 <- plot_waffle(dfsub, "Ethnicity") +
  scale_fill_manual(values = c("Latino" = "#FF5722", "Asian" = "#FF9800", "Caucasian" = "#2196F3", "African" = "#4CAF50", "Jewish" = "#9C27B0", "Mixed" = "#795548"))

p8 <- plot_waffle(dfsub, "Screen_Resolution", title = "Screen Resolution") +
  scale_fill_pizza_d() +
  guides(fill = "none")

p9 <- plot_waffle(dfsub, "Device_OS", title = "Device OS") +
  scale_fill_bluebrown_d()

# p10 <- plot_waffle(dfsub, "Screen_Refresh") +
#   scale_fill_viridis_d()


(p1 / p2 / p3) | (p4 / p5 / p9) | (p6 / p7 / p8)

Models

Delboeuf

Model Selection

data <- filter(df, Illusion_Type == "Delboeuf")

# mm <- model.matrix(Error ~ abs(Illusion_Strength) * Illusion_Effect, data=data)
# head(as.data.frame(mm))

best_models(data)
##                        Name   BIC R2_marginal       BF     Side  Model
## 1  DIFF_sqrt--STRENGTH_sqrt 11023      0.4205          6.62e-05 Errors
## 2  DIFF_cbrt--STRENGTH_sqrt 11026      0.4125  = 0.195 6.63e-05 Errors
## 3   DIFF_log--STRENGTH_sqrt 11038      0.4335  < 0.001 6.73e-05 Errors
## 4  DIFF_sqrt--STRENGTH_cbrt 11046      0.4345  < 0.001 6.91e-05 Errors
## 5   DIFF_sqrt--STRENGTH_log 11048      0.4116  < 0.001 6.62e-05 Errors
## 6   DIFF_sqrt--STRENGTH_log 11418      0.0710          2.70e-01     RT
## 7    DIFF_log--STRENGTH_log 11421      0.0709  = 0.175 2.66e-01     RT
## 8   DIFF_cbrt--STRENGTH_log 11422      0.0708  = 0.122 2.72e-01     RT
## 9  DIFF_sqrt--STRENGTH_sqrt 11432      0.0703  < 0.001 2.60e-01     RT
## 10      DIFF_sqrt--STRENGTH 11433      0.0702  < 0.001 2.72e-01     RT

Error Rate

Descriptive

plot_descriptive_err(data)

Model Specification

formula <- brms::bf(
  Error ~ sqrt(Illusion_Difference) * sqrt(abs(Illusion_Strength)) * Illusion_Effect +
    (1 + sqrt(Illusion_Difference) * sqrt(abs(Illusion_Strength)) * Illusion_Effect | Participant),
  family = "bernoulli"
)

# brms::get_prior(formula, data = data)
# brms::validate_prior(formula)

model_delboeuf_err <- brms::brm(formula,
  data = data,
  refresh = 0
)

# parameters::parameters(model_delboeuf_err)

save(model_delboeuf_err, file="models/model_delboeuf_err.Rdata")

load("models/model_delboeuf_err.RData")

Model Visualization

p_delboeuf_err <- plot_model_err(data, model_delboeuf_err)
p_delboeuf_err

Reaction Time

Descriptive

data <- filter(df, Illusion_Type == "Delboeuf", Error == 0)

plot_descriptive_rt(data)

Model Specification

formula <- brms::bf(
  RT ~ sqrt(Illusion_Difference) * sqrt(abs(Illusion_Strength)) * Illusion_Effect +
    (1 + sqrt(Illusion_Difference) * sqrt(abs(Illusion_Strength)) * Illusion_Effect | Participant),
  sigma ~ 1 + (1 | Participant),
  beta ~ 1 + (1 | Participant),
  family = "exgaussian"
)


# brms::get_prior(formula, data = data)
# plot(seq(-5, 5, length.out=100), brms::dstudent_t(seq(-5, 5, length.out=100), df = 3, 0.1, 0.5))
priors <- c(
  # Intercept
  prior("student_t(3, 0, 0.1)", class = "Intercept"),
  prior("student_t(3, 0.35, 0.1)", class = "sd", coef="Intercept", group="Participant"),
  # Default parameters
  prior("student_t(3, 0, 3)", class = "b"),
  prior("student_t(3, 0.35, 0.5)", class = "sd", group="Participant")
  
  # prior("student_t(3, 0, 1)", class="Intercept", dpar = "sigma"),
  # prior("student_t(3, 0.1, 1)", class="sd", dpar = "sigma")
  
  # prior("student_t(3, 0.1, 0.5)", class = "ndt")
  ) |> 
  brms::validate_prior(formula, data = data)

# bayestestR::model_to_priors(model_delboeuf_rt, scale_multiply=3)

model_delboeuf_rt <- brms::brm(formula,
  data = data,
  refresh = 50,
  init=0
)

save(model_delboeuf_rt, file="models/model_delboeuf_rt.Rdata")

load("models/model_delboeuf_rt.RData")

plot_ppcheck(model_delboeuf_rt)
# performance::check_model(model_delboeuf_rt)
# parameters::parameters(model)

Model Visualization

p_delboeuf_rt <- plot_model_rt(data, model_delboeuf_rt)
p_delboeuf_rt

Visualization

p_delboeuf <- plot_all(data, p_delboeuf_err, p_delboeuf_rt)
p_delboeuf

Ebbinghaus

Model Selection

data <- filter(df, Illusion_Type == "Ebbinghaus")

best_models(data)
##                       Name  BIC R2_marginal       BF     Side  Model
## 1  DIFF_sqrt--STRENGTH_log 8184      0.5750          8.73e-09 Errors
## 2   DIFF_log--STRENGTH_log 8185      0.5948  = 0.454 9.21e-09 Errors
## 3      DIFF_sqrt--STRENGTH 8193      0.5638  = 0.010 8.77e-09 Errors
## 4       DIFF_log--STRENGTH 8198      0.5790  = 0.001 9.33e-09 Errors
## 5  DIFF_cbrt--STRENGTH_log 8200      0.5655  < 0.001 8.83e-09 Errors
## 6      DIFF_sqrt--STRENGTH 8336      0.0943          2.31e-01     RT
## 7  DIFF_sqrt--STRENGTH_log 8338      0.0942  = 0.447 2.33e-01     RT
## 8      DIFF_cbrt--STRENGTH 8339      0.0942  = 0.320 2.32e-01     RT
## 9  DIFF_cbrt--STRENGTH_log 8340      0.0941  = 0.154 2.34e-01     RT
## 10      DIFF_log--STRENGTH 8345      0.0938  = 0.013 2.32e-01     RT

Error Rate

Descriptive

plot_descriptive_err(data)

Model Specification

formula <- brms::bf(
  Error ~ sqrt(Illusion_Difference) * logmod(abs(Illusion_Strength)) * Illusion_Effect +
    (1 + sqrt(Illusion_Difference) * logmod(abs(Illusion_Strength)) * Illusion_Effect | Participant),
  family = "bernoulli"
)

# brms::get_prior(formula, data = data)
# brms::validate_prior(formula)

model_ebbinghaus_err <- brms::brm(formula,
  data = data,
  refresh = 0
)

# parameters::parameters(model_delboeuf_err)

save(model_ebbinghaus_err, file="models/model_ebbinghaus_err.Rdata")

load("models/model_ebbinghaus_err.RData")

Model Visualization

p_ebbinghaus_err <- plot_model_err(data, model_ebbinghaus_err)
p_ebbinghaus_err

Reaction Time

Descriptive

data <- filter(df, Illusion_Type == "Ebbinghaus", Error == 0)

plot_descriptive_rt(data)

Model Specification

formula <- brms::bf(
  RT ~ sqrt(Illusion_Difference) * logmod(abs(Illusion_Strength)) * Illusion_Effect +
    (1 + sqrt(Illusion_Difference) * logmod(abs(Illusion_Strength)) * Illusion_Effect | Participant),
  sigma ~ 1 + (1 | Participant),
  beta ~ 1 + (1 | Participant),
  family = "exgaussian"
)


model_ebbinghaus_rt <- brms::brm(formula,
  data = data,
  refresh = 0,
  init = 0
)

save(model_ebbinghaus_rt, file="models/model_ebbinghaus_rt.Rdata")

load("models/model_ebbinghaus_rt.RData")

plot_ppcheck(model_ebbinghaus_rt)
# performance::check_model(model_delboeuf_rt)
# parameters::parameters(model)

Model Visualization

p_ebbinghaus_rt <- plot_model_rt(data, model_ebbinghaus_rt)
p_ebbinghaus_rt

Rod and Frame

Model Selection

data <- filter(df, Illusion_Type == "Rod-Frame")

best_models(data)
##                        Name   BIC R2_marginal       BF     Side  Model
## 1    DIFF_log--STRENGTH_log 10354      0.4875          0.000146 Errors
## 2   DIFF_sqrt--STRENGTH_log 10357      0.4893  = 0.208 0.000148 Errors
## 3   DIFF_log--STRENGTH_sqrt 10367      0.4845  = 0.002 0.000146 Errors
## 4  DIFF_sqrt--STRENGTH_sqrt 10370      0.4865  < 0.001 0.000147 Errors
## 5   DIFF_log--STRENGTH_cbrt 10374      0.4905  < 0.001 0.000149 Errors
## 6    DIFF_log--STRENGTH_log  8688      0.0900          0.076184     RT
## 7   DIFF_log--STRENGTH_sqrt  8689      0.0899  = 0.435 0.076600     RT
## 8        DIFF_log--STRENGTH  8695      0.0895  = 0.028 0.079896     RT
## 9   DIFF_log--STRENGTH_cbrt  8700      0.0891  = 0.002 0.077029     RT
## 10  DIFF_sqrt--STRENGTH_log  8701      0.0891  = 0.001 0.077129     RT

Error Rate

Descriptive

plot_descriptive_err(data)

Model Specification

formula <- brms::bf(
  Error ~ logmod(Illusion_Difference) * logmod(abs(Illusion_Strength)) * Illusion_Effect +
    (1 + logmod(Illusion_Difference) * logmod(abs(Illusion_Strength)) * Illusion_Effect | Participant),
  family = "bernoulli"
)

# brms::get_prior(formula, data = data)
# brms::validate_prior(formula)

model_rodframe_err <- brms::brm(formula,
  data = data,
  refresh = 0
)
## Running MCMC with 4 parallel chains...
## 
## Chain 3 finished in 1585.6 seconds.
## Chain 1 finished in 1610.1 seconds.
## Chain 2 finished in 1616.1 seconds.
## Chain 4 finished in 1618.5 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 1607.6 seconds.
## Total execution time: 1618.8 seconds.

# parameters::parameters(model_delboeuf_err)

save(model_rodframe_err, file="models/model_rodframe_err.Rdata")

load("models/model_rodframe_err.RData")

Model Visualization

p_rodframe_err <- plot_model_err(data, model_rodframe_err)
p_rodframe_err

Reaction Time

Descriptive

data <- filter(df, Illusion_Type == "Rod-Frame", Error == 0)

plot_descriptive_rt(data)

Model Specification

formula <- brms::bf(
  RT ~ logmod(Illusion_Difference) * logmod(abs(Illusion_Strength)) * Illusion_Effect +
    (1 | Participant),
  sigma ~ 1 + (1 | Participant),
  beta ~ 1 + (1 | Participant),
  family = "exgaussian"
)


model_rodframe_rt <- brms::brm(formula,
  data = data,
  refresh = 0,
  init = 0
)

save(model_rodframe_rt, file="models/model_rodframe_rt.Rdata")

load("models/model_rodframe_rt.RData")

plot_ppcheck(model_rodframe_rt)
# performance::check_model(model_delboeuf_rt)
# parameters::parameters(model)

Model Visualization

p_rodframe_rt <- plot_model_rt(data, model_rodframe_rt)
p_rodframe_rt

Vertical-Horizontal

Model Selection

data <- filter(df, Illusion_Type == "Vertical-Horizontal")

best_models(data)
##                        Name   BIC R2_marginal       BF  Side  Model
## 1  DIFF_sqrt--STRENGTH_sqrt 10047       0.579          0.128 Errors
## 2  DIFF_cbrt--STRENGTH_sqrt 10054       0.578  = 0.027 0.129 Errors
## 3   DIFF_log--STRENGTH_sqrt 10061       0.581  < 0.001 0.129 Errors
## 4       DIFF--STRENGTH_sqrt 10071       0.582  < 0.001 0.129 Errors
## 5  DIFF_sqrt--STRENGTH_cbrt 10143       0.586  < 0.001 0.129 Errors
## 6   DIFF_log--STRENGTH_sqrt  7234       0.128          0.875     RT
## 7  DIFF_sqrt--STRENGTH_sqrt  7234       0.128  = 0.956 0.888     RT
## 8       DIFF--STRENGTH_sqrt  7237       0.128  = 0.189 0.871     RT
## 9  DIFF_cbrt--STRENGTH_sqrt  7238       0.128  = 0.089 0.893     RT
## 10  DIFF_log--STRENGTH_cbrt  7247       0.127  = 0.001 0.880     RT

Error Rate

Descriptive

plot_descriptive_err(data)

Model Specification

formula <- brms::bf(
  Error ~ sqrt(Illusion_Difference) * sqrt(abs(Illusion_Strength)) * Illusion_Effect +
    (1 + sqrt(Illusion_Difference) * sqrt(abs(Illusion_Strength)) * Illusion_Effect | Participant),
  family = "bernoulli"
)

# brms::get_prior(formula, data = data)
# brms::validate_prior(formula)

model_verticalhorizontal_err <- brms::brm(formula,
  data = data,
  refresh = 0
)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 finished in 1850.4 seconds.
## Chain 4 finished in 1869.3 seconds.
## Chain 2 finished in 1916.9 seconds.
## Chain 3 finished in 1929.4 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 1891.5 seconds.
## Total execution time: 1929.7 seconds.

# parameters::parameters(model_verticalhorizontal_err)

save(model_verticalhorizontal_err, file="models/model_verticalhorizontal_err.Rdata")

load("models/model_verticalhorizontal_err.RData")

Model Visualization

p_verticalhorizontal_err <- plot_model_err(data, model_verticalhorizontal_err)
p_verticalhorizontal_err

Reaction Time

Descriptive

data <- filter(df, Illusion_Type == "Vertical-Horizontal", Error == 0)

plot_descriptive_rt(data)

Model Specification

formula <- brms::bf(
  RT ~ sqrt(Illusion_Difference) * sqrt(abs(Illusion_Strength)) * Illusion_Effect +
    (1 | Participant),
  # sigma ~ t2(Illusion_Difference * Illusion_Strength) +
  #   (1 | Participant),
  family = "exgaussian"
)


model_verticalhorizontal_rt <- brms::brm(formula,
  data = data,
  refresh = 0
)

save(model_verticalhorizontal_rt, file="models/model_verticalhorizontal_rt.Rdata")

load("models/model_verticalhorizontal_rt.RData")

plot_ppcheck(model_verticalhorizontal_rt)
# performance::check_model(model_verticalhorizontal_rt)
# parameters::parameters(model)

Model Visualization

p_verticalhorizontal_rt <- plot_model_rt(data, model_verticalhorizontal_rt)
p_verticalhorizontal_rt

Zöllner

Model Selection

data <- filter(df, Illusion_Type == "Zöllner")

best_models(data)
##                        Name   BIC R2_marginal       BF    Side  Model
## 1       DIFF_cbrt--STRENGTH  9367      0.4069          0.71416 Errors
## 2        DIFF_log--STRENGTH  9375      0.4154  = 0.013 0.71320 Errors
## 3       DIFF_sqrt--STRENGTH  9376      0.4170  = 0.009 0.71262 Errors
## 4            DIFF--STRENGTH  9521      0.4468  < 0.001 0.70663 Errors
## 5  DIFF_cbrt--STRENGTH_sqrt  9772      0.3878  < 0.001 0.70402 Errors
## 6       DIFF_cbrt--STRENGTH 13383      0.0690          0.00455     RT
## 7        DIFF_log--STRENGTH 13384      0.0689  = 0.545 0.00467     RT
## 8       DIFF_sqrt--STRENGTH 13396      0.0683  = 0.001 0.00483     RT
## 9  DIFF_cbrt--STRENGTH_sqrt 13420      0.0670  < 0.001 0.00442     RT
## 10  DIFF_log--STRENGTH_sqrt 13420      0.0669  < 0.001 0.00455     RT

Error Rate

Descriptive

plot_descriptive_err(data)

Model Specification

formula <- brms::bf(
  Error ~ cbrtmod(Illusion_Difference) * abs(Illusion_Strength) * Illusion_Effect +
    (1 + cbrtmod(Illusion_Difference) * abs(Illusion_Strength) * Illusion_Effect | Participant),
  family = "bernoulli"
)

# brms::get_prior(formula, data = data)
# brms::validate_prior(formula)

model_zollner_err <- brms::brm(formula,
  data = data,
  refresh = 0
)
## Running MCMC with 4 parallel chains...
## 
## Chain 3 finished in 1926.9 seconds.
## Chain 1 finished in 1964.4 seconds.
## Chain 2 finished in 1973.6 seconds.
## Chain 4 finished in 1980.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 1961.3 seconds.
## Total execution time: 1980.3 seconds.

# parameters::parameters(model_zollner_err)

save(model_zollner_err, file="models/model_zollner_err.Rdata")

load("models/model_zollner_err.RData")

Model Visualization

p_zollner_err <- plot_model_err(data, model_zollner_err)
p_zollner_err

Reaction Time

Descriptive

data <- filter(df, Illusion_Type == "Zöllner", Error == 0)

plot_descriptive_rt(data)

Model Specification

formula <- brms::bf(
  RT ~ cbrtmod(Illusion_Difference) * abs(Illusion_Strength) * Illusion_Effect +
    (1 | Participant),
  # sigma ~ t2(Illusion_Difference * Illusion_Strength) +
  #   (1 | Participant),
  family = "exgaussian"
)


model_zollner_rt <- brms::brm(formula,
  data = data,
  refresh = 0
)

save(model_zollner_rt, file="models/model_zollner_rt.Rdata")

load("models/model_zollner_rt.RData")

plot_ppcheck(model_zollner_rt)
# performance::check_model(model_zollner_rt)
# parameters::parameters(model)

Model Visualization

p_zollner_rt <- plot_model_rt(data, model_zollner_rt)
p_zollner_rt

White

data <- filter(df, Illusion_Type == "White")

best_models(data)
##                        Name  BIC R2_marginal       BF    Side  Model
## 1  DIFF_sqrt--STRENGTH_sqrt 5821      0.7677          0.00175 Errors
## 2       DIFF--STRENGTH_sqrt 5827      0.7809  = 0.045 0.00165 Errors
## 3  DIFF_cbrt--STRENGTH_sqrt 5834      0.7644  = 0.002 0.00182 Errors
## 4   DIFF_log--STRENGTH_sqrt 5861      0.7615  < 0.001 0.00193 Errors
## 5   DIFF_sqrt--STRENGTH_log 5862      0.7970  < 0.001 0.00200 Errors
## 6        DIFF_log--STRENGTH 9590      0.0903          0.01015     RT
## 7       DIFF_cbrt--STRENGTH 9600      0.0897  = 0.008 0.01027     RT
## 8       DIFF_sqrt--STRENGTH 9610      0.0892  < 0.001 0.01035     RT
## 9   DIFF_log--STRENGTH_sqrt 9644      0.0866  < 0.001 0.01129     RT
## 10           DIFF--STRENGTH 9648      0.0869  < 0.001 0.01048     RT

Error Rate

Descriptive

plot_descriptive_err(data)

Model Specification

formula <- brms::bf(
  Error ~ Illusion_Difference * sqrt(abs(Illusion_Strength)) * Illusion_Effect +
    (1 + Illusion_Difference * sqrt(abs(Illusion_Strength)) * Illusion_Effect | Participant),
  family = "bernoulli"
)

# brms::get_prior(formula, data = data)
# brms::validate_prior(formula)

model_white_err <- brms::brm(formula,
  data = data,
  refresh = 0
)
## Start sampling
## Running MCMC with 4 parallel chains...
## 
## Chain 1 finished in 2677.9 seconds.
## Chain 4 finished in 2991.8 seconds.
## Chain 3 finished in 3058.0 seconds.
## Chain 2 finished in 3117.4 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 2961.3 seconds.
## Total execution time: 3117.8 seconds.

# parameters::parameters(model_white_err)

save(model_white_err, file="models/model_white_err.Rdata")

load("models/model_white_err.RData")

Model Visualization

p_white_err <- plot_model_err(data, model_white_err)
p_white_err

Reaction Time

Descriptive

data <- filter(df, Illusion_Type == "White", Error == 0)

plot_descriptive_rt(data)

Model Specification

formula <- brms::bf(
  RT ~ Illusion_Difference * sqrt(abs(Illusion_Strength)) * Illusion_Effect +
    (1 | Participant),
  # sigma ~ t2(Illusion_Difference * Illusion_Strength) +
  #   (1 | Participant),
  family = "exgaussian"
)


model_white_rt <- brms::brm(formula,
  data = data,
  refresh = 0
)

save(model_white_rt, file="models/model_white_rt.Rdata")

load("models/model_white_rt.RData")

plot_ppcheck(model_white_rt)
# performance::check_model(model_white_rt)
# parameters::parameters(model)

Model Visualization

p_white_rt <- plot_model_rt(data, model_white_rt)
p_white_rt

Müller-Lyer

Model Selection

data <- filter(df, Illusion_Type == "Müller-Lyer")

best_models(data)
##                        Name  BIC R2_marginal       BF     Side  Model
## 1  DIFF_sqrt--STRENGTH_sqrt 7273       0.732          5.25e-17 Errors
## 2   DIFF_log--STRENGTH_sqrt 7273       0.730  = 0.738 4.28e-17 Errors
## 3       DIFF--STRENGTH_sqrt 7291       0.730  < 0.001 4.32e-17 Errors
## 4  DIFF_cbrt--STRENGTH_sqrt 7292       0.734  < 0.001 6.54e-17 Errors
## 5    DIFF_log--STRENGTH_log 7299       0.752  < 0.001 5.76e-17 Errors
## 6       DIFF_sqrt--STRENGTH 9593       0.154          1.64e-09     RT
## 7       DIFF_cbrt--STRENGTH 9593       0.154  = 0.953 1.62e-09     RT
## 8  DIFF_cbrt--STRENGTH_sqrt 9602       0.153  = 0.007 6.02e-10     RT
## 9  DIFF_sqrt--STRENGTH_sqrt 9604       0.153  = 0.004 6.19e-10     RT
## 10       DIFF_log--STRENGTH 9605       0.153  = 0.002 1.86e-09     RT

Error Rate

Descriptive

plot_descriptive_err(data, side = "updown")

Model Specification

formula <- brms::bf(
  Error ~ sqrt(Illusion_Difference) * sqrt(abs(Illusion_Strength)) * Illusion_Effect +
    (1 + sqrt(Illusion_Difference) * sqrt(abs(Illusion_Strength)) * Illusion_Effect | Participant),
  family = "bernoulli"
)

# brms::get_prior(formula, data = data)
# brms::validate_prior(formula)

model_mullerlyer_err <- brms::brm(formula,
  data = data,
  refresh = 0
)
## Running MCMC with 4 parallel chains...
## 
## Chain 3 finished in 2054.6 seconds.
## Chain 1 finished in 2112.0 seconds.
## Chain 4 finished in 2506.1 seconds.
## Chain 2 finished in 2610.9 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 2320.9 seconds.
## Total execution time: 2611.3 seconds.

# parameters::parameters(model_white_err)

save(model_mullerlyer_err, file="models/model_mullerlyer_err.Rdata")

load("models/model_mullerlyer_err.RData")

Model Visualization

p_mullerlyer_err <- plot_model_err(data, model_mullerlyer_err)
p_mullerlyer_err

Reaction Time

Descriptive

data <- filter(df, Illusion_Type == "Müller-Lyer", Error == 0)

plot_descriptive_rt(data, side = "updown")

Model Specification

formula <- brms::bf(
  RT ~ sqrt(Illusion_Difference) * sqrt(abs(Illusion_Strength)) * Illusion_Effect +
    (1 | Participant),
  # sigma ~ t2(Illusion_Difference * Illusion_Strength) +
  #   (1 | Participant),
  family = "exgaussian"
)


model_mullerlyer_rt <- brms::brm(formula,
  data = data,
  refresh = 0
)

save(model_mullerlyer_rt, file="models/model_mullerlyer_rt.Rdata")
# parameters::parameters(model_white_err)

load("models/model_mullerlyer_rt.RData")

plot_ppcheck(model_mullerlyer_rt)
# performance::check_model(model_mullerlyer_rt)
# parameters::parameters(model)

Model Visualization

p_mullerlyer_rt <- plot_model_rt(data, model_mullerlyer_rt)
p_mullerlyer_rt

Ponzo

Model Selection

data <- filter(df, Illusion_Type == "Ponzo")

best_models(data)
##                        Name   BIC R2_marginal       BF     Side  Model
## 1       DIFF_cbrt--STRENGTH  8250       0.563          0.735137 Errors
## 2       DIFF_sqrt--STRENGTH  8254       0.570  = 0.135 0.737259 Errors
## 3        DIFF_log--STRENGTH  8304       0.590  < 0.001 0.742920 Errors
## 4            DIFF--STRENGTH  8339       0.599  < 0.001 0.745206 Errors
## 5  DIFF_cbrt--STRENGTH_sqrt  8368       0.586  < 0.001 0.761796 Errors
## 6       DIFF_cbrt--STRENGTH 10267       0.118          0.000664     RT
## 7       DIFF_sqrt--STRENGTH 10279       0.117  = 0.003 0.000631     RT
## 8        DIFF_log--STRENGTH 10325       0.115  < 0.001 0.000579     RT
## 9            DIFF--STRENGTH 10357       0.114  < 0.001 0.000566     RT
## 10 DIFF_cbrt--STRENGTH_sqrt 10362       0.113  < 0.001 0.001393     RT

Error Rate

Descriptive

plot_descriptive_err(data, side = "updown")

Model Specification

formula <- brms::bf(
  Error ~ cbrtmod(Illusion_Difference) * abs(Illusion_Strength) * Illusion_Effect +
    (1 + cbrtmod(Illusion_Difference) * abs(Illusion_Strength) * Illusion_Effect | Participant),
  family = "bernoulli"
)

# brms::get_prior(formula, data = data)
# brms::validate_prior(formula)

model_ponzo_err <- brms::brm(formula,
  data = data,
  refresh = 0
)
## Running MCMC with 4 parallel chains...
## 
## Chain 2 finished in 2044.8 seconds.
## Chain 4 finished in 2077.4 seconds.
## Chain 3 finished in 2093.1 seconds.
## Chain 1 finished in 2096.9 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 2078.1 seconds.
## Total execution time: 2097.2 seconds.

# parameters::parameters(model_ponzo_err)

save(model_ponzo_err, file="models/model_ponzo_err.Rdata")

load("models/model_ponzo_err.RData")

Model Visualization

p_ponzo_err <- plot_model_err(data, model_ponzo_err)
p_ponzo_err

Reaction Time

Descriptive

data <- filter(df, Illusion_Type == "Ponzo", Error == 0)

plot_descriptive_rt(data, side = "updown")

Model Specification

formula <- brms::bf(
  RT ~ cbrtmod(Illusion_Difference) * abs(Illusion_Strength) * Illusion_Effect +
    (1 | Participant),
  # sigma ~ t2(Illusion_Difference * Illusion_Strength) +
  #   (1 | Participant),
  family = "exgaussian"
)


model_ponzo_rt <- brms::brm(formula,
  data = data,
  refresh = 0
)

save(model_ponzo_rt, file="models/model_ponzo_rt.Rdata")
# parameters::parameters(model_white_err)

load("models/model_ponzo_rt.RData")

plot_ppcheck(model_ponzo_rt)
# performance::check_model(model_ponzo_rt)
# parameters::parameters(model)

Model Visualization

p_ponzo_rt <- plot_model_rt(data, model_ponzo_rt)
p_ponzo_rt

Poggendorff

Model Selection

data <- filter(df, Illusion_Type == "Poggendorff")

best_models(data)
##                        Name  BIC R2_marginal       BF   Side  Model
## 1  DIFF_cbrt--STRENGTH_sqrt 8259      0.5190          0.6545 Errors
## 2  DIFF_sqrt--STRENGTH_sqrt 8281      0.5193  < 0.001 0.6607 Errors
## 3       DIFF_cbrt--STRENGTH 8325      0.4986  < 0.001 0.6401 Errors
## 4       DIFF_sqrt--STRENGTH 8348      0.4973  < 0.001 0.6463 Errors
## 5   DIFF_log--STRENGTH_sqrt 8380      0.5184  < 0.001 0.6746 Errors
## 6       DIFF_cbrt--STRENGTH 8461      0.0778          0.0251     RT
## 7       DIFF_sqrt--STRENGTH 8480      0.0766  < 0.001 0.0246     RT
## 8  DIFF_cbrt--STRENGTH_sqrt 8520      0.0740  < 0.001 0.0268     RT
## 9        DIFF_log--STRENGTH 8534      0.0734  < 0.001 0.0233     RT
## 10 DIFF_sqrt--STRENGTH_sqrt 8540      0.0728  < 0.001 0.0263     RT

Error Rate

Descriptive

plot_descriptive_err(data, side = "updown")

Model Specification

formula <- brms::bf(
  Error ~ cbrtmod(Illusion_Difference) * sqrt(abs(Illusion_Strength)) * Illusion_Effect +
    (1 + cbrtmod(Illusion_Difference) * sqrt(abs(Illusion_Strength)) * Illusion_Effect | Participant),
  family = "bernoulli"
)

# brms::get_prior(formula, data = data)
# brms::validate_prior(formula)

model_poggendorff_err <- brms::brm(formula,
  data = data,
  refresh = 0
)
## Running MCMC with 4 parallel chains...
## 
## Chain 4 finished in 3086.8 seconds.
## Chain 3 finished in 3150.5 seconds.
## Chain 2 finished in 3172.7 seconds.
## Chain 1 finished in 3187.4 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 3149.4 seconds.
## Total execution time: 3187.6 seconds.

# parameters::parameters(model_poggendorff_err)

save(model_poggendorff_err, file="models/model_poggendorff_err.Rdata")

load("models/model_poggendorff_err.RData")

Model Visualization

p_poggendorff_err <- plot_model_err(data, model_poggendorff_err)
p_poggendorff_err

Reaction Time

Descriptive

data <- filter(df, Illusion_Type == "Poggendorff", Error == 0)

plot_descriptive_rt(data, side = "updown")

Model Specification

formula <- brms::bf(
  RT ~ cbrtmod(Illusion_Difference) * sqrt(abs(Illusion_Strength)) * Illusion_Effect +
    (1 | Participant),
  # sigma ~ 1 + (1 | Participant),
  family = "exgaussian"
)


model_poggendorff_rt <- brms::brm(formula,
  data = data,
  refresh = 0
)

save(model_ponzo_rt, file="models/model_poggendorff_rt.Rdata")

load("models/model_poggendorff_rt.RData")

plot_ppcheck(model_poggendorff_rt)
# performance::check_model(model_poggendorff_rt)
# parameters::parameters(model)

Model Visualization

p_poggendorff_rt <- plot_model_rt(data, model_poggendorff_rt)
p_poggendorff_rt

Contrast

Model Selection

data <- filter(df, Illusion_Type == "Contrast")

best_models(data)
##                        Name   BIC R2_marginal       BF     Side  Model
## 1  DIFF_sqrt--STRENGTH_sqrt  7837       0.610          4.70e-19 Errors
## 2       DIFF--STRENGTH_sqrt  7842       0.621  = 0.084 4.55e-19 Errors
## 3  DIFF_cbrt--STRENGTH_sqrt  7847       0.607  = 0.008 5.12e-19 Errors
## 4   DIFF_sqrt--STRENGTH_log  7859       0.658  < 0.001 8.41e-19 Errors
## 5        DIFF--STRENGTH_log  7864       0.675  < 0.001 8.08e-19 Errors
## 6  DIFF_sqrt--STRENGTH_sqrt 10454       0.122          2.42e-02     RT
## 7  DIFF_cbrt--STRENGTH_sqrt 10454       0.122  = 0.824 2.43e-02     RT
## 8   DIFF_log--STRENGTH_sqrt 10457       0.122  = 0.211 2.45e-02     RT
## 9       DIFF--STRENGTH_sqrt 10462       0.122  = 0.013 2.44e-02     RT
## 10      DIFF_sqrt--STRENGTH 10479       0.121  < 0.001 1.94e-02     RT

Error Rate

Descriptive

data <- filter(df, Illusion_Type == "Contrast")

plot_descriptive_err(data, side = "updown")

Model Specification

formula <- brms::bf(
  Error ~ sqrt(Illusion_Difference) * sqrt(abs(Illusion_Strength)) * Illusion_Effect +
    (1 + sqrt(Illusion_Difference) * sqrt(abs(Illusion_Strength)) * Illusion_Effect | Participant),
  family = "bernoulli"
)

# brms::get_prior(formula, data = data)
# brms::validate_prior(formula)

model_contrast_err <- brms::brm(formula,
  data = data,
  refresh = 0
)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 finished in 3483.4 seconds.
## Chain 2 finished in 3506.1 seconds.
## Chain 4 finished in 3551.5 seconds.
## Chain 3 finished in 3555.2 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 3524.0 seconds.
## Total execution time: 3555.6 seconds.

# parameters::parameters(model_contrast_err)

save(model_contrast_err, file="models/model_contrast_err.Rdata")

load("models/model_contrast_err.RData")

Model Visualization

p_contrast_err <- plot_model_err(data, model_contrast_err)
p_contrast_err

Reaction Time

Descriptive

data <- filter(df, Illusion_Type == "Contrast", Error == 0)

plot_descriptive_rt(data, side = "updown")

Model Specification

formula <- brms::bf(
  RT ~ sqrt(Illusion_Difference) * sqrt(abs(Illusion_Strength)) * Illusion_Effect +
    (1 | Participant),
  # sigma ~ 1 + (1 | Participant),
  family = "exgaussian"
)


model_contrast_rt <- brms::brm(formula,
  data = data,
  refresh = 0
)

save(model_contrast_rt, file="models/model_contrast_rt.Rdata")

load("models/model_contrast_rt.RData")

plot_ppcheck(model_contrast_rt)
# performance::check_model(model_contrast_rt)
# parameters::parameters(model)

Model Visualization

p_contrast_rt <- plot_model_rt(data, model_contrast_rt)
p_contrast_rt

Individual Scores

Reliability

random <- rbind(
  extract_random(model_delboeuf_err, "Delboeuf"),
  # extract_random(model_delboeuf_rt, "Delboeuf")
  extract_random(model_ebbinghaus_err, "Ebbinghaus"),
  # extract_random(model_ebbinghaus_rt, "Ebbinghaus"),
  extract_random(model_rodframe_err, "Rod-Frame"),
  # extract_random(model_rodframe_rt, "Rod-Frame"),
  extract_random(model_verticalhorizontal_err, "Vertical-Horizontal"),
  # extract_random(model_verticalhorizontal_rt, "Vertical-Horizontal"),
  extract_random(model_zollner_err, "Zöllner"),
  # extract_random(model_zollner_rt, "Zöllner"),
  extract_random(model_white_err, "White"),
  # extract_random(model_white_rt, "White"),
  extract_random(model_mullerlyer_err, "Müller-Lyer"),
  # extract_random(model_mullerlyer_rt, "Müller-Lyer"),
  extract_random(model_ponzo_err, "Ponzo"),
  # extract_random(model_ponzo_rt, "Ponzo"),
  extract_random(model_poggendorff_err, "Poggendorff"),
  # extract_random(model_poggendorff_rt, "Poggendorff"),
  extract_random(model_contrast_err, "Contrast")
  # extract_random(model_contrast_rt, "Contrast")
) |> 
  filter(
    # !str_detect(Parameter, "Intercept_"),
    !str_detect(Parameter, "Cong_"),
    !str_detect(Parameter, "Null_"))


random |>
  # group_by(Illusion_Type, Parameter) |> 
  # standardize(select="Value") |> 
  # ungroup() |> 
  # mutate(order = mean(Value)) |> 
  # ungroup() |> 
  # mutate(Participant = fct_reorder(Participant, order)) |> 
  ggplot(aes(y = Participant, x = Value)) +
  # ggdist::stat_slabinterval(aes(fill = Illusion_Type, color=Illusion_Type), alpha=0.5) +
  ggdist::stat_pointinterval(aes(fill = Illusion_Type, color=Illusion_Type), alpha=0.5) +
  # coord_cartesian(xlim=c(-4, 4)) +
  facet_wrap(~Parameter, scales = "free_x") +
  theme(axis.text.y = element_blank())



random <- random |>
  mutate(Parameter = paste0(Illusion_Type, "_", Parameter),
         Parameter = format_illusionName(Parameter)) |>
  group_by(Participant, Parameter) |>
  summarize(Score = median(Value), .groups = "drop") |>
  ungroup() |>
  group_by(Parameter) |> 
  mutate(Score = Score / sd(Score)) |>  # Scale
  ungroup()  |> 
  pivot_wider(names_from = Parameter, values_from = Score)

# plot(estimate_density(select(random, -Participant)))

Relationship with empirical descriptors

# Average when incongruent
empirical <- df |>
  filter(Illusion_Effect == "Incongruent") |> 
  group_by(Participant, Illusion_Type) |>
  summarize(
    # n = n(),
    Error_General = sum(Error, na.rm=TRUE) / n(),
    RT_Mean_General = mean(RT, na.rm=TRUE),
    RT_SD_General = sd(RT, na.rm=TRUE),
  ) |>
  ungroup() |> 
  mutate(Illusion_Type = format_illusionName(Illusion_Type))  |> 
  pivot_wider(names_from = c("Illusion_Type"), values_from = c("Error_General", "RT_Mean_General", "RT_SD_General")) |> 
  arrange(Participant)

# Average for 3 higher level of illusion strength
empirical <- df |>
  filter(Illusion_Effect == "Incongruent") |> 
  group_by(Participant, Illusion_Type, Illusion_Strength) |>
  summarize(
    # n = n(),
    Error_Strong = sum(Error, na.rm=TRUE) / n(),
    RT_Mean_Strong = mean(RT, na.rm=TRUE),
    RT_SD_Strong = sd(RT, na.rm=TRUE),
  ) |>
  slice((n()-2):n()) |> 
  group_by(Participant, Illusion_Type) |>
  summarize(
    # n = n(),
    Error_Strong = mean(Error_Strong, na.rm=TRUE),
    RT_Mean_Strong = mean(RT_Mean_Strong, na.rm=TRUE),
    RT_SD_Strong = mean(RT_SD_Strong, na.rm=TRUE),
  ) |> 
  ungroup() |> 
  mutate(Illusion_Type = format_illusionName(Illusion_Type)) |> 
  pivot_wider(names_from = c("Illusion_Type"), values_from = c("Error_Strong", "RT_Mean_Strong", "RT_SD_Strong")) |> 
  arrange(Participant) |> 
  cbind(empirical)

# Correlation
r <- data.frame()
for(ill in c("Delboeuf", "Ebbinghaus", 
             # "RodFrame", 
             "VerticalHorizontal", 
             # "Zollner", 
             # "White", 
             "MullerLyer", 
             "Ponzo", 
             # "Poggendorff", 
             "Contrast")) {
  r <- correlation::correlation(
    data = random |> 
      arrange(Participant) |> 
      data_select(select=ill, regex=TRUE),
    data2 = empirical |>
      data_select(select=ill, regex=TRUE) 
    ) |> 
    mutate(Illusion_Type = ill,
           Parameter1 = clean_name(Parameter1),
           Parameter2 = clean_name(Parameter2)) |> 
    select(Illusion_Type, Parameter1, Parameter2, r, p) |> 
    rbind(r)
}

r |>
  ggplot(aes(x = Parameter2, y = Parameter1)) +
  geom_tile(aes(fill = r, alpha = p)) +
  scale_alpha_continuous(range = c(1, 0)) +
  scale_fill_gradient2(low = "#2196F3", mid = "white", high = "#F44336", midpoint = 0, limit = c(-1, 1), space = "Lab", name = "Correlation", guide = "legend") +
  facet_wrap(~Illusion_Type, scales = "free")

Structure

Correlation

cor <- correlation::correlation(random, redundant = TRUE)

p_data <- cor |>
  mutate(
    Parameter1 = clean_name(Parameter1),
    Parameter2 = clean_name(Parameter2)
  ) |>
  correlation::cor_sort(hclust_method = "ward.D2") |>
  cor_lower() |>
  mutate(
    Text = insight::format_value(r, zap_small = TRUE, digits = 3),
    Text = str_replace(str_remove(Text, "^0+"), "^-0+", "-"),
    Parameter2 = fct_rev(Parameter2)
  )


p_data |>
  ggplot(aes(x = Parameter2, y = Parameter1)) +
  geom_tile(aes(fill = r)) +
  # geom_text(aes(label = Text), size = 2) +
  scale_fill_gradient2(low = "#2196F3", mid = "white", high = "#F44336", midpoint = 0, limit = c(-1, 1), space = "Lab", name = "Correlation", guide = "legend") +
  scale_x_discrete(expand = c(0, 0)) +
  scale_y_discrete(expand = c(0, 0)) +
  labs(title = "Correlation Matrix", x = NULL, y = NULL) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 90, hjust = 1),
    plot.title = element_text(hjust = 0.5, face = "bold"),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
  )

Factor Analysis

First Level

data <- select(random, -Participant)

rez <- parameters::n_factors(data)
plot(rez)


fa <- parameters::factor_analysis(data,
  cor = as.matrix(cor),
  n = 8,
  rotation = "oblimin",
  fm = "mle",
  sort = TRUE
)

# psych::omega(data, nfactors=3)
# fa <- psych::fa.multi(as.matrix(cor), nfactors = 3, nfact2 = 1, n.obs = nrow(random))
# psych::fa.multi.diagram(fa)

Second Level

data <- predict(fa)

rez <- parameters::n_factors(data, rotation = "oblimin")
plot(rez)

fa2 <- parameters::factor_analysis(data,
  n = 2,
  rotation = "varimax",
  fm = "mle",
  sort = TRUE
)

# psych::schmid(attributes(fa)$model, nfactors = 1)

Visualization

library(ggraph)

dat <- rbind(
  attributes(fa)$loadings_long,
  attributes(fa2)$loadings_long |>
    mutate(Component = str_replace(Component, "ML", "G"))
) |>
  tidygraph::as_tbl_graph() |>
  tidygraph::activate("nodes")


# Get order of indices
# idx <- correlation::cor_sort(cor, hclust_method = "ward.D2") |>
#   summary() |>
#   pull(Parameter)
# idx <- sort(fa$Variable)
idx <- fa$Variable

get_x <- function(name, x) {
  case_when(
    name == "G1" ~ 0/1,
    name == "G2" ~ 1/1,
    name == "ML1" ~ 0/4,
    name == "ML2" ~ 1/4,
    name == "ML3" ~ 2/4,
    name == "ML4" ~ 3/4,
    name == "ML5" ~ 4/4,
    TRUE ~ x)
}

dat |>
  mutate(
    x = sapply(as.data.frame(dat)$name, function(x) {
      if (x %in% idx) {
        (which(idx == x) - 1) / (length(idx) - 1)
      } else {
        NA
      }
    })
  ) |>
  mutate(
    x = get_x(name, x),
    y = case_when(
    str_detect(name, "G") ~ 2,
    str_detect(name, "ML") ~ 1,
    TRUE ~ 0),
    angle = ifelse(name %in% idx, 90, 0),
    hjust = ifelse(name %in% idx, 1.01, 0),
    name = clean_name(name),
    Illusion_Type = str_split_fixed(name, " - ", n=2)[, 1]) |> 
  tidygraph::activate("edges") |> 
  filter(abs(Loading) > 0.2) |> 
  ggraph(layout = "nicely") +
  geom_edge_link(aes(edge_width = Loading, edge_colour = as.factor(sign(Loading)))) +
  geom_node_point(aes(color = Illusion_Type), size=3) +
  geom_node_text(aes(label = name, angle = angle, hjust=hjust)) +
  scale_y_continuous(expand = expansion(c(.6, .1))) +
  scale_edge_width_continuous(range=c(0.1, 1)) +
  scale_edge_color_manual(values=c("#F44336", "#2196F3")) +
  scale_color_manual(values=c("Delboeuf" = "#F44336", Ebbinghaus = "#E91E63",
                              "Rod-Frame" = "#9C27B0")) +
  guides(edge_colour = "none") +
  ggraph::theme_graph()

SEM

library(lavaan)
library(tidySEM)


names <- names(select(random, -Participant))
data <- setNames(as.data.frame(str_split_fixed(names, "_", 3)), c("Illusion_Type", "Parameter", "Model"))
data$Name <- names

# Model 1
lvl1 <- paste(data$Illusion_Type, "=~", data$Name)
lvl2 <- paste("G =~", unique(data$Illusion_Type))
resid <- as.data.frame(t(utils::combn(x=unique(data$Illusion_Type), 2)))
resid <- paste(resid$V1, "~~", resid$V2)[1:2]



m1 <- sem(model = paste0(c(lvl1, lvl2, resid), collapse = "\n"),
           data  = random)


# Model 2
lvl1 <- paste(data$Model, "=~", data$Name)
lvl2 <- paste("G =~", unique(data$Model))

m2 <- sem(model = paste0(c(lvl1, lvl2), collapse = "\n"),
           data  = random)

# Model 3
lvl1 <- paste(data$Parameter, "=~", data$Name)
lvl2 <- paste("G =~", unique(data$Parameter))

m3 <- sem(model = paste0(c(lvl1, lvl2), collapse = "\n"),
           data  = random)

anova(m1, m2, m3)

summary(m1, standardize = TRUE)

tidySEM::graph_sem(model = m1)

# as_tbl_graph(fit, standardize = TRUE)

# Save
dfsub <- dfsub |> 
  merge(random, by = "Participant") |> 
  cbind(data, predict(fa2, names=c("G1", "G2")))

write.csv("../data/study3.csv")

Illusion Game Validation (Study 2)

2022-08-17

Introduction

Methods

Exclusions

Error Rate

Reaction Time Distribution

Error Rate per Illusion Block

Reaction Time per Trial

Participants

Models

Delboeuf

Model Selection

Error Rate

Descriptive

Model Specification

Model Visualization

Reaction Time

Descriptive

Model Specification

Model Visualization

Visualization

Ebbinghaus

Model Selection

Error Rate

Descriptive

Model Specification

Model Visualization

Reaction Time

Descriptive

Model Specification

Model Visualization

Rod and Frame

Model Selection

Error Rate

Descriptive

Model Specification

Model Visualization

Reaction Time

Descriptive

Model Specification

Model Visualization

Vertical-Horizontal

Model Selection

Error Rate

Descriptive

Model Specification

Model Visualization

Reaction Time

Descriptive

Model Specification

Model Visualization

Zöllner

Model Selection

Error Rate

Descriptive

Model Specification

Model Visualization

Reaction Time

Descriptive

Model Specification

Model Visualization

White

Error Rate

Descriptive

Model Specification

Model Visualization

Reaction Time

Descriptive

Model Specification

Model Visualization

Müller-Lyer

Model Selection

Error Rate

Descriptive

Model Specification

Model Visualization

Reaction Time

Descriptive

Model Specification